﻿/*******************************************************************
Copyright (c) 2011 강성진 이재범

*  Permission is hereby granted, free of charge, to any person
*  obtaining a copy of this software and associated documentation
*  files (the "Software"), to deal in the Software without
*  restriction, including without limitation the rights to use,
*  copy, modify, merge, publish, distribute, sublicense, and/or sell
*  copies of the Software, and to permit persons to whom the
*  Software is furnished to do so, subject to the following
*  conditions:
* 
*  The above copyright notice and this permission notice shall be
*  included in all copies or substantial portions of the Software.
* 
*  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
*  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
*  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
*  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
*  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
*  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
*  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
*  OTHER DEALINGS IN THE SOFTWARE.
*******************************************************************/

#include "common.h"
#include<cmath>
using namespace std;

static double eps = 10e-9;

void outer(double ** arr, double a[],double b[], int n)
	//arr에 a*b해주기
{
	for(int q = 0; q < n; q++)
	{
		for(int p = 0; p < n; p++)
		{
			arr[p][q] = a[q]*b[p];
		}
	}
}
double** made_matrix(int n)
	//n크기의 행렬생성
{
	double **arr;
	arr = new double* [n];
	for(int i = 0;i < n;i++)
	{
		arr[i] = new double[n];
	}
	return arr;
}
int find_pivot(double **p, int a)
	//p행렬(a*a)의 피봇 기준 (절대값이 가장 큰 것을 찾는다)
{
	double i = p[0][0];
	int row_max = 0;
	for(int b = 0; b < a; b++)
	{
		if(abs(i) < abs(p[0][b]))
		{
			row_max = b;
			i = p[0][b];
		}
	}
	return row_max;
}
void permutation(double**p, int a, int b,int n, int& plmn)
	// permutation matrix p 에서 row a switch with row b at n*n
{
	if(a == b)
		return;

	double *arr;

	arr=new double[n];

	for(int i=0;i<n;i++)
		arr[i] = p[i][a];

	for(int i=0;i<n;i++)
		p[i][a] = p[i][b];

	for(int i=0;i<n;i++)
		p[i][b] = arr[i];

	delete [] arr;
	plmn *= -1;
}
int calculate_determinant(int** relation_matrix, int dim)
{
	int plmn = 1;
	double **mat;
	double *pivot;

	pivot= new double[dim];

	int probe = 0;

	while(probe < dim)
	{
		// 행렬 입력 받는 부분
		if(probe == 0)
		{
			mat = new double* [dim];
			for(int i = 0;i < dim;i++)
			{
				mat[i] = new double [dim];
			}

			for(int q = 0; q < dim; q++)
			{
				for(int p = 0; p < dim; p++)
					mat[p][q] = relation_matrix[p][q];
			}
		}

		permutation(mat, 0, find_pivot(mat, dim - probe), dim - probe, plmn);

		pivot[probe]=mat[0][0];

		if(abs(pivot[probe]) <= eps)
			return 0;

		double *v;//v 받을 배열 생성
		v = new double[dim - probe - 1];

		double *w;//w 받을 배열 생성
		w = new double[dim - probe - 1];

		for(int i = 1;i < dim - probe; i++)//v,w 초기화
		{
			v[i-1] = mat[0][i];
			w[i-1] = mat[i][0];
		}

		probe += 1;

		double ** vwT;// vw 외적 저장공간
		vwT = new double*[dim-probe];

		for(int i = 0;i < dim - probe;i++)
		{
			vwT[i] = new double[dim - probe];
		}

		outer(vwT, v, w, dim - probe);//외적 저장

		double ** a_prime; //a'
		a_prime = new double* [dim - probe];

		for(int i = 0;i < dim - probe;i++)
		{
			a_prime[i] = new double [dim - probe];
		}

		for(int j = 0; j < dim - probe;j++)
		{
			for(int i = 0;i < dim - probe;i++)
			{
				a_prime[i][j] = mat[i+1][j+1];
			}
		}

		for(int i = 0;i < dim - probe;i++)
		{
			delete mat[i];
		}

		delete mat;
		mat = new double* [dim - probe];//mat에 (A'-vwT)/aii대입
		for(int i = 0;i < dim - probe;i++)
		{
			mat[i] = new double [dim - probe];
		}

		for(int j = 0;j < dim - probe;j++)
		{
			for(int i = 0;i < dim - probe;i++)
			{
				mat[i][j] = a_prime[i][j] - (vwT[i][j]/pivot[probe - 1]);
			}

		}
		if(probe == (dim - 1))
		{
			pivot[dim - 1] = mat[0][0];

			if(abs(pivot[dim - 1]) <= eps)
				return 0;
			break;
		}

		delete v;
		delete w;

		for(int i = 0;i < dim - probe;i++)
		{
			delete vwT[i];
		}

		delete vwT;
	}

	double determinant = 1;

	for(int i = 0; i < dim;i++)
		determinant *= pivot[i];

	return static_cast<int>(determinant * plmn);
}

